#!/usr/bin/env python
import numpy as np
import itk
import scipy.ndimage as ndi
import sys
import os
import pylab as pb

def vas_seg(fname, low_th, high_th, start, end, outputfile):
#Get segment the lung vasculature and grab the largest num of components
	reader = itk.ImageFileReader.ISS3.New()
        reader.SetFileName(fname)
	btf = itk.BinaryThresholdImageFilter.ISS3ISS3.New()
	btf.SetInput(reader.GetOutput())
	btf.SetLowerThreshold(low_th)
	btf.SetUpperThreshold(high_th)
	ccf = itk.ConnectedComponentImageFilter.ISS3IUL3.New()
	ccf.SetInput(btf.GetOutput())
	rcf = itk.RelabelComponentImageFilter.IUL3IUS3.New()
	rcf.SetInput(ccf.GetOutput())
        rcf.Update()
	tf = itk.ThresholdImageFilter.IUS3.New()
	tf.SetInput(rcf.GetOutput())
	tf.SetLower(start)
	tf.SetUpper(end)
	tf.Update()
        writer= itk.ImageFileWriter.IUS3.New()
        writer.SetFileName(outputfile)
        writer.SetInput(tf.GetOutput())
        writer.Update()
        return rcf.GetSizeOfObjectsInPixels()
def center(fname):
       	reader = itk.ImageFileReader.IUC3.New()
        reader.SetFileName(fname)
        img = reader.GetOutput()
        array = itk.PyBuffer.IUC3.GetArrayFromImage(img)
        idx = np.where(array>0)
        #print idx[0]
        #print "mean value is", np.mean(idx[0]), np.mean(idx[1]), np.mean(idx[2]), 
        #print "var is ", np.std(idx[0]), np.std(idx[1]), np.std(idx[2])
        #stat = [np.mean(idx[0]), np.mean(idx[1]), np.mean(idx[2]), np.std(idx[0]), np.std(idx[1]), np.std(idx[2])] 
        stat = [np.std(idx[0]), np.std(idx[1]), np.std(idx[2])] 
        return stat 
def main():
       if len(sys.argv) < 6 :
            print " Usage : vas_seg fname, low_th, high_th, start, end"
            sys.exit(1)
       low_th = int(sys.argv[2])
       high_th = int(sys.argv[3]) 
       tmp = os.path.splitext(sys.argv[1])
       outputfile = tmp[0] +"__" + str(low_th) + "__"+str(high_th) + tmp[1] 
       size = vas_seg(sys.argv[1], low_th, high_th, int(sys.argv[4]), int(sys.argv[5]), outputfile)
       print size[0:5]
       pb.plot(size[0:5])
       pb.show() 
if __name__ == '__main__':
    main()

